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Abstract. Restoration of images degraded by spatially varying blurs is an issue of increasing importance 
in the context of photography, satellite or microscopy imaging. One of the main difficulty to solve this 
problem comes from the huge dimensions of the blur matrix. It prevents the use of naive approaches for 
performing matrix-vector multiplications. In this paper, we propose to approximate the blur operator by 
a matrix sparse in the wavelet domain. We justify this approach from a mathematical point of view and 
investigate the approximation quality numerically. We finish by showing that the sparsity pattern of the 
matrix can be pre-defined, which is central in tasks such as blind deconvolution. 



1. Introduction 

This problem of image restoration in the presence of spatially varying blurs drew a considerable attention 
in the context of satellite imaging with the Hubble space telescope [IJ. It is now becoming increasingly 
important in various imaging modalities such as 3D fluorescence microscopy 13. In this paper, we aim at 
providing new computational tools to tackle this problem. We consider a blurring operator H inQ. cW^, 
defined for any u G L^(t2) as the following integral operator, 

Vx G Hu{x) = / K{x,y)u{y)dy, (1) 

where f2 C is the image domain. The function ^ is a spatially varying kernel defining the Point 
Spread Function (PSF) at each location x ^Q.. 

The naive computation of a product Hu is simply obtained by direct discretization of ([T]) and costs 
O {N^^^ arithmetic operations for an image of dimension d with N pixels in each dimension. It is 
unsuitable for large scale problems encountered in imaging. Previous works attempted to reduce this 
computational cost by approximating H using either tensor products or piecewise convolutions EITIIH. 
Unfortunately, these strategies do not have a property that is highly desirable in numerical analysis: the 
ability to approximate the original operator with an arbitrary precision. In practice, we observed that 
they are often too coarse for complex situations encountered in microscopy imaging. 

In a recent work (H, we proposed to approximate H by an operator diagonal in the wavelet domain. 
However, this approximation was shown to be too coarse for some practical applications. Our main 
contribution in this paper consists of showing that the diagonal approximation can be replaced by a sparse 
approximation and lead to arbitrarily close approximations. More precisely, we show that \\H — ^0^* || , 
can be made arbitrarily small if ^* denotes a forward wavelet transform and 5 is a sparse matrix that 



contains only 0{N^) non zero coefficients instead of N^^. As a practical consequence, matrix- vector 
products using the proposed approach costs O (d\og{N)N^) arithmetic operations, which is doable even 
for very large images. We also show that the sparsity pattern of S can be known a priori, which is a 
highly desirable property in the case of blind deconvolution, where the operator should be inferred from 
the blurred images. 

In section [2} we recall the theoretical results motivating this approximation. Section |3.1[ contains 
preliminary restoration results in order to illustrate theoretical results and to identify a number of 
coefficients offering a compromise between the storage costs and the restoration quality. In this section, 
the kept coefficients are obtained by simple thresholding and can have arbitrary coordinates. In section 
|3.2[ we finally show that the sparsity pattern of the matrix can be pre-defined and that this approach 
leads to almost equivalent results to a thresholding. 

2. Theoretical motivations 

Sparse approximations of integral operators have been theoretically analyzed in ||6l|7l. Surprisingly 
this approach was never applied to the approximation of spatially varying blur operators. The closest 
application found in the literature deals with foveation [8 |, but it is not adapted to a large class of kernels. 
We define an orthonormal wavelet basis of L^(R) as the family of functions: 

Each y^fj^yn is a dilated and translated of the mother wavelet \|/, 

Let us recall a typical result that motivate the proposed approach. We stick to the one-dimensional 
case for the ease of exposition. Since // is a linear operator in a Hilbert space, it can be written as 
H — ^0^*, where : ^ is the matrix representation of the blur operator in the wavelet domain. 
Matrix is characterized by the coefficients: 

Theorem 1 (Decay of Qj^m,k,n - lUl). Supposing that: 

• the wavelets are compactly supported with supp(\\fj^rn) — ^j^m, 

• the wavelets have M vanishing moments, 

• the operator H belongs to the class of Calderon-Zygmund operators, meaning that K satisfies: 

{ \K{x,y)\<-^ 

I l-^~3^l 
3Cm > 0, such that ^x^y G fl, < ^ (2) 

\\d^Kix,y)\ + \d^Kix,y)\<- m<1 



then the coefficients Qj^m,k,n satisfy the following estimate: 

|e, .„| < C„mi„ ( M M I ' ( . (3) 

[14,^1 |^j,m|J \ dlSt{Ij^rnJk,n) J 

with Cm a constant depending on M and the wavelets, dist{Ij^mih,n) denotes the distance between the 
two supports of the wavelets and y^k.n- 

Moreover, if the kernel K is compactly supported, then for sufficiently large \m — n\. 



^j,m,k,n — 0. 



(4) 



This results ensures that the coefficients away from the diagonal fastly decay to zero. A more precise 
analysis shows that O (^dN^ log{N)) non zero coefficients allow to approximate H with an arbitrary 
precision |6|. Figure [T] shows the matrix associated to a typical spatially varying blur operator in 
log scale. It is readily seen that the decrease of the coefficients away from the diagonal is extremely fast. 
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Figure 1. An example of a matrix in a 
logjo scale for a 256 x 256 images. The matrix 
contains about 4.3 billion coefficients. 



Figure 2. The point spread function of the 
blur used to generate the matrix on the left. 
The kernel is a Gaussian with equal variances 
increasing in the vertical direction. 



In the next section, we reconstruct degraded images from the knowledge of various thresholded 
matrices. 



3. Restoration results 

3.1. Restoration using thresholded matrices 

In a first experiment, the matrix is thresholded in order to keep a number of coefficients being a multiple 
of (the number of pixels of the image). That is, we define the thresholded matrix 07^^ by zeroing all but 
the Tk = kN^ largest coefficients of 0. This experiment allows the identification of an admissible number 
of coefficients offering a proper compromise between the storage costs and the image restoration quality. 
To achieve this experiment, a full matrix is numerically computed and then successively thresholded 
for various ^ G N*. Thereafter, the image is restored assuming the following classical degradation model: 



where v is the degraded image observed (see Figure Bq), u is the image to restore (see Figure 3a), and H 



is a blurring operator corresponding to the PSF displayed in Figure [3d| A standard TV-L2 optimization 
problem is solved to restore the image u: 



Find w* G argmin TV{u), (5) 

mgR""^ , 1 1 ^0r^^* w- V 1 1 2 <c^n 

where TV{u) is the isotropic total variation of u. 

Figure|4]displays restored images for various values of k and the sparsity pattern of 0^^. As expected, 
considering more coefficients improves the deblurring results. However, in the presence of noise, the 
two images restored using the full matrix and the thresholded matrix with k = 20, provide very similar 
results. Therefore, it is possible to store only a small amount of coefficients and obtain results close to 
the ones obtained having a full (unrealistic) knowledge of the blur matrix. Furthermore, these results are 
obtained drastically faster, see the speed-up Figure |4j 



Remark: it is near impossible to solve the deblurring problem with an exact operator due to the 
numerical burden of computing matrix-vector products in large dimensions. We thus chose a spatially 
varying blur that is sparse in the wavelet domain in order to be able to make comparisons between an 
exact computation and an approximated one. Speed up times might thus be much higher for other kinds 
of blurs. 




(a) Original Image (b) Blurred Image (c) Degraded Image a = (d) PSF 

2.10-2, 57V/? = 16.54dB. 



Figure 3. Images involved in the deblurring problem 




(a) Non-zeros coefficient for ^ = 1 




(d) Deblurred image for = 1 . SNR = 
17.12dB, Speed-up = 22. 



(b) Non-zeros coefficient for k = 20 




(e) Deblurred image for k = 20. 
SNR = 18.69dB, Speed-up = 17. 




(f) Deblurred image with the full 0, 
kc:^ 1550, SNR =n.93dB. 



Figure 4. Deblurred images using thresholded matrices 



3.2. Restoration using pre-defined sparsity patterns 

We are currently developping strategies to solve blind deconvolution problems using the proposed 
approximation. In this setting, it is not possible to threshold the exact (unknown) matrix to obtain 
a sparsity pattern. The knowledge of such a structure is primordial to approximate K using prior 
information. Our aim in this section is to show that sparsity patterns can be pre-defined, allowing 
convincing restoration results. 

Our main observation is that for a given wavelet only its neighbours in the same scale, the 
finer scale and the coarser one lead to significant correlation coefficients {H\\fj^fn^^k,n)' This concept of 
neighbours is also supported by Theorem[T] Thus, sparsity patterns can be defined using notions of intra 
and inter- scale neighbourhoods. 

For each node of the wavelet decomposition quadtree, a neighbourhood is defined. It describes which 
correlation coefficients shall be preserved to generate the sparsity pattern. Figure |5] illustrates this idea. 
In this example, we consider a wavelet transform of depth 2. A neighbourhood 9^ is associated to each 
sub-band. For instance 0\(x can be represented as follows: 



/2 2 2 2 2 2 2 2\ 

/ / / / I h V d 

-1 1 

\o -1 1 oy 



scale 
oriention 

vertical translation 
horizontal translation 



where / stands for the low frequency wavelet, /z, v, J for the horizontal, vertical and diagonal orientations 
respectively. Figure [6] illustrates this neighbourhood fA^ for a given wavelet at the center of the image. 




Figure 5. Quadtree corresponding to a wavelet 

transform of depth 2 ^'S"'"^ ^' Illustration of the neighbourhood 

IAq. In black, the given wavelet at the center 
of the image, and in red all the preserved 
neighbouring wavelets. 



In order to validate this notion of neighbourhood, we considered two different scenarii. The first one 

/ same scale \ 



consists of setting for all /, lA^ -■ 



all orientations 





The second one consists of setting for all /, 



( same scale same scale same scale same scale same scale \ 

all orientations all orientations all orientations all orientations all orientations 
0-1100 

\ -1 1 / 



Figure[7]displays restored images using sparsity patterns generated from these two neighbourhoods. This 
experiment shows that sparsity patterns defined from the correlations between neighbouring wavelets are 
relevant for our aim: approximating and interpolating the blurring operator. 
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(a) Non-zero coefficients of pattern 1 




(c) Deblurred image with pattern 1, 
kc:^3,SNR=16.14dB. 
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(b) Non-zero coefficients of pattern 2 




(d) Deblurred image with pattern 2, 
kc:^l5,SNR=lS30dB. 



Figure 7. Deblurred images using matrices defined by sparsity patterns 



Conclusion 

In this paper we investigated the use of sparse matrices in the wavelet domain in order to approximate 
spatially varying blur operators. We observed that very large compression factors still allow to obtain 
near optimal reconstruction results. Compared to previously proposed approach, this technique allows 
to approximate the exact operator with an arbitrary precision. It is thus suited to noiseless and noisy 
problems. We also showed that pre-established sparsity patterns can be used to efficiently approximate 
the blurring integral operator. This property is central in order to tackle blind deconvolution problems. 
We are currently working on this topic. 
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